####Analyses with the full sample


######GAM-Models (Primary Model)#######

#Run all models for the full country sample in one loop
model_GAM_Fullcountry = datCompleteNew %>% 
  filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5")) %>% 
  filter(!(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty" | dv_name == "DemImportant")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = bam(dv_value ~ age_grp1 + age_grp2 + waves.f + s(born_adult) + sex, data = .)
  )

#Save the regression output for the Shiny App
source("11_RegTables.R")


name_flag <- model_GAM_Fullcountry$dv_name
country_flag <- model_GAM_Fullcountry$country

#Create new data for predictions
datBind <- list()

for(i in 1:length(name_flag)){
  datBind[[i]] <- datCompleteNew %>% 
    filter(waves %in% c("EVS1","EVS2","EVS3","EVS4","EVS5")) %>% 
    filter(!(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty" | dv_name == "DemImportant")) %>% 
    filter(dv_name == name_flag[i]) %>% 
    filter(country == country_flag[i]) %>% 
    drop_na(age_grp1, age_grp2, waves.f, born_adult, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, waves.f, born_adult, age_grp1, age_grp2, sex) %>% 
    droplevels() %>% 
    dplyr::select(waves.f, born_adult, age_grp1, age_grp2, sex)
}


predList <- list()

for(i in 1:length(name_flag)){
  predList[[i]] <- predict(model_GAM_Fullcountry$fitCountry[[i]], datBind[[i]],  type = "link", se.fit = TRUE)
}

preds <- list()

for(i in 1:length(name_flag)){
  preds[[i]] <- transform(cbind(data.frame(predList[[i]]), datBind[[i]]))
}


#Calculate critical value for simultaneous intervals
rmvn <- function(n, mu, sig) { ## MVN random deviates
  L <- mroot(sig)
  m <- ncol(L)
  t(mu + L %*% matrix(rnorm(m*n), m, n))
}


crit <- list()

for(i in 1:length(name_flag)){
  
  Vb <- vcov(model_GAM_Fullcountry$fitCountry[[i]])
  
  pred <- predict(model_GAM_Fullcountry$fitCountry[[i]], se.fit = TRUE)
  
  se.fit <- pred$se.fit
  
  set.seed(42)
  N <- 1000
  
  BUdiff <- rmvn(N, mu = rep(0, nrow(Vb)), sig = Vb)
  
  Cg <- predict(model_GAM_Fullcountry$fitCountry[[i]], type = "lpmatrix")
  simDev <- Cg %*% t(BUdiff)
  
  absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))
  
  masd <- apply(absDev, 2L, max)
  
  crit[[i]] <- quantile(masd, prob = 0.95, type = 8)
}



#Aggregate predictions
predBornAdult <- list()

for(i in 1:length(name_flag)){
  
  predBornAdult[[i]] <- preds[[i]] %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}


predPeriod <- list()

for(i in 1:length(name_flag)){
  predPeriod[[i]] <- preds[[i]] %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}

names(predBornAdult) <- c(paste(name_flag, as.character(country_flag), sep = ";"))

predBornAdult <- bind_rows(predBornAdult, .id = "dv_name")

predBornAdult <- predBornAdult %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")


names(predPeriod) <- c(paste(name_flag, as.character(country_flag), sep = ";"))

predPeriod <- bind_rows(predPeriod, .id = "dv_name")

predPeriod <- predPeriod %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "EVS1" = 1981, "EVS2" = 1990, "EVS3" = 1999, "EVS4" = 2008, "EVS5" = 2017))

#####Plotting#####


#Create plot of predictions for all countries and all DV and store them as png
pg_list <- list()

for(i in 1:length(predList)){
  
  message(i)
  
  pg_list[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult %>% filter(dv_name == name_flag[i], country == country_flag[i])
    pp_period.data <- predPeriod %>% filter(dv_name == name_flag[i], country == country_flag[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(-1, 1, 0.4), limits = c(-1,1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      geom_errorbar(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, ymin = LL, ymax = UL), colour = cols[1], width=.5, position=position_dodge(0.1)) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020), limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}

for(i in 1:24){
  pg_list[[i]] <- pg_list[[i]] + scale_y_continuous(name="", limits = c(-1.1,1.3)) 
}

for(i in 25:length(predList)){
  pg_list[[i]] <- pg_list[[i]] + scale_y_continuous(name="", breaks=seq(0, 1, 0.2), limits = c(-0.1,1.1)) 
}


for(i in seq(1,length(predList), by = 24)){
  agg_png(file=paste0(wdShinyPlots, "plotsFull_GAM_cohort_",i,".png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list[[i]], pg_list[[i+1]], pg_list[[i+2]], pg_list[[i+3]], pg_list[[i+4]], pg_list[[i+5]], pg_list[[i+6]], pg_list[[i+7]], pg_list[[i+8]], pg_list[[i+9]], pg_list[[i+10]], pg_list[[i+11]], pg_list[[i+12]], pg_list[[i+13]], pg_list[[i+14]], pg_list[[i+15]], pg_list[[i+16]], pg_list[[i+17]], pg_list[[i+18]], pg_list[[i+19]], pg_list[[i+20]], pg_list[[i+21]], pg_list[[i+22]], pg_list[[i+23]], nrow = 4))
  
  dev.off()
}



#########Subsample ###############

#Run the GAM models for the smaller country sample (Democracy Important & Essential variables)
model_GAM_Subcountry_DemImportant = datCompleteNew %>% 
  filter(country %in% countrySub_DemImportant) %>% 
  filter((dv_name == "DemImportant")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = bam(dv_value ~ age_grp1 + age_grp2 + waves.f + s(born_adult) + sex, data = .)
  )

model_GAM_Subcountry_Essentials = datCompleteNew %>% 
  filter(country %in% countrySub_Essentials) %>% 
  filter((dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty")) %>% 
  group_by(dv_name, country) %>%
  do(fitCountry = bam(dv_value ~ age_grp1 + age_grp2 + waves.f + s(born_adult) + sex, data = .)
  )

#Save regression tables for Shiny App
source("13_RegTables_Subsample.R")



name_flag_DemImportant <- model_GAM_Subcountry_DemImportant$dv_name
country_flag_DemImportant <- model_GAM_Subcountry_DemImportant$country

name_flag_Essentials <- model_GAM_Subcountry_Essentials$dv_name
country_flag_Essentials <- model_GAM_Subcountry_Essentials$country


#Create new data for predictions
datBind_DemImportant <- list()
datBind_Essentials <- list()


for(i in 1:length(name_flag_DemImportant)){
  datBind_DemImportant[[i]] <- datCompleteNew %>% 
    filter(waves %in% c("EVS5", "WVS5", "WVS6")) %>% 
    filter(dv_name == "DemImportant") %>% 
    filter(dv_name == name_flag_DemImportant[i]) %>% 
    filter(country == country_flag_DemImportant[i]) %>% 
    drop_na(age_grp1, age_grp2, waves.f, born_adult, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, waves.f, born_adult, age_grp1, age_grp2, sex) %>% 
    droplevels() %>% 
    dplyr::select(waves.f, born_adult, age_grp1, age_grp2, sex)
}

for(i in 1:length(name_flag_Essentials)){
  datBind_Essentials[[i]]  <- datCompleteNew %>% 
    filter(waves %in% c("EVS5", "WVS5", "WVS6")) %>% 
    filter(dv_name == "democracy_freeElection" | dv_name == "democracy_ArmyTakeOver" | dv_name == "democracy_protectLiberty") %>% 
    filter(dv_name == name_flag_Essentials[i]) %>% 
    filter(country == country_flag_Essentials[i]) %>% 
    drop_na(age_grp1, age_grp2, waves.f, born_adult, dv_value, sex) %>% 
    dplyr::select(dv_name, dv_value, waves.f, born_adult, age_grp1, age_grp2, sex) %>% 
    droplevels() %>% 
    dplyr::select(waves.f, born_adult, age_grp1, age_grp2, sex)
}

predList_DemImportant <- list()
predList_Essentials <- list()


for(i in 1:length(name_flag_DemImportant)){
  predList_DemImportant[[i]] <- predict(model_GAM_Subcountry_DemImportant$fitCountry[[i]], datBind_DemImportant[[i]],  type = "link", se.fit = TRUE)
}

for(i in 1:length(name_flag_Essentials)){
  predList_Essentials[[i]] <- predict(model_GAM_Subcountry_Essentials$fitCountry[[i]], datBind_Essentials[[i]],  type = "link", se.fit = TRUE)
}

preds_DemImportant <- list()
preds_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  preds_DemImportant[[i]] <- transform(cbind(data.frame(predList_DemImportant[[i]]), datBind_DemImportant[[i]]))
}

for(i in 1:length(name_flag_Essentials)){
  preds_Essentials[[i]] <- transform(cbind(data.frame(predList_Essentials[[i]]), datBind_Essentials[[i]]))
}


#Calculate critical values for simultaneous intervals
rmvn <- function(n, mu, sig) { ## MVN random deviates
  L <- mroot(sig)
  m <- ncol(L)
  t(mu + L %*% matrix(rnorm(m*n), m, n))
}


crit_DemImportant <- list()
crit_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  
  Vb <- vcov(model_GAM_Subcountry_DemImportant$fitCountry[[i]])
  
  pred <- predict(model_GAM_Subcountry_DemImportant$fitCountry[[i]], se.fit = TRUE)
  
  se.fit <- pred$se.fit
  
  set.seed(42)
  N <- 1000
  
  BUdiff <- rmvn(N, mu = rep(0, nrow(Vb)), sig = Vb)
  
  Cg <- predict(model_GAM_Subcountry_DemImportant$fitCountry[[i]], type = "lpmatrix")
  simDev <- Cg %*% t(BUdiff)
  
  absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))
  
  masd <- apply(absDev, 2L, max)
  
  crit_DemImportant[[i]] <- quantile(masd, prob = 0.95, type = 8)
}

for(i in 1:length(name_flag_Essentials)){
  
  Vb <- vcov(model_GAM_Subcountry_Essentials$fitCountry[[i]])
  
  pred <- predict(model_GAM_Subcountry_Essentials$fitCountry[[i]], se.fit = TRUE)
  
  se.fit <- pred$se.fit
  
  set.seed(42)
  N <- 1000
  
  BUdiff <- rmvn(N, mu = rep(0, nrow(Vb)), sig = Vb)
  
  Cg <- predict(model_GAM_Subcountry_Essentials$fitCountry[[i]], type = "lpmatrix")
  simDev <- Cg %*% t(BUdiff)
  
  absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))
  
  masd <- apply(absDev, 2L, max)
  
  crit_Essentials[[i]] <- quantile(masd, prob = 0.95, type = 8)
}




#Aggregate predictions
predBornAdult_DemImportant <- list()
predBornAdult_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  
  predBornAdult_DemImportant[[i]] <- preds_DemImportant[[i]] %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}


for(i in 1:length(name_flag_Essentials)){
  
  predBornAdult_Essentials[[i]] <- preds_Essentials[[i]] %>% 
    group_by(born_adult) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(born_adult, fit, UL, LL)
  
}


predPeriod_DemImportant <- list()
predPeriod_Essentials <- list()

for(i in 1:length(name_flag_DemImportant)){
  predPeriod_DemImportant[[i]] <- preds_DemImportant[[i]] %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}

for(i in 1:length(name_flag_Essentials)){
  predPeriod_Essentials[[i]] <- preds_Essentials[[i]] %>% 
    group_by(waves.f) %>% 
    mutate(m = n(),
           mean_fit = mean(fit),
           S_sq = (fit - mean_fit)^2/(m - 1)) %>%
    replace_na(list(S_sq = 0)) %>% 
    summarize(fit = mean(fit),
              S_sq = mean(S_sq),
              se.fit = mean(se.fit^2) + S_sq,
              se.fit = sqrt(se.fit)) %>% 
    ungroup() %>% 
    mutate(UL = fit + crit[[i]] * se.fit,
           LL = fit - crit[[i]] * se.fit) %>% 
    dplyr::select(waves.f, fit, UL, LL)
  
}


names(predBornAdult_DemImportant) <- c(paste(name_flag_DemImportant, as.character(country_flag_DemImportant), sep = ";"))
names(predBornAdult_Essentials) <- c(paste(name_flag_Essentials, as.character(country_flag_Essentials), sep = ";"))

predBornAdult_DemImportant <- bind_rows(predBornAdult_DemImportant, .id = "dv_name")
predBornAdult_Essentials <- bind_rows(predBornAdult_Essentials, .id = "dv_name")

predBornAdult_DemImportant <- predBornAdult_DemImportant %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")
predBornAdult_Essentials <- predBornAdult_Essentials %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";")


names(predPeriod_DemImportant) <- c(paste(name_flag_DemImportant, as.character(country_flag_DemImportant), sep = ";"))
names(predPeriod_Essentials) <- c(paste(name_flag_Essentials, as.character(country_flag_Essentials), sep = ";"))

predPeriod_DemImportant <- bind_rows(predPeriod_DemImportant, .id = "dv_name")
predPeriod_Essentials <- bind_rows(predPeriod_Essentials, .id = "dv_name")

predPeriod_DemImportant <- predPeriod_DemImportant %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "WVS5" = 2007, "WVS6" = 2012, "EVS5" = 2017))

predPeriod_Essentials<- predPeriod_Essentials %>% 
  tidyr::separate(dv_name, c("dv_name", "country"), ";") %>% 
  mutate(year = recode(waves.f, "WVS5" = 2007, "WVS6" = 2012, "EVS5" = 2017))


#####Plotting######

pg_list_DemImportant <- list()
pg_list_Essentials <- list()

for(i in 1:length(predList_DemImportant)){
  
  message(i)
  
  pg_list_DemImportant[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult_DemImportant %>% filter(dv_name == name_flag_DemImportant[i], country == country_flag_DemImportant[i])
    pp_period.data <- predPeriod_DemImportant %>% filter(dv_name == name_flag_DemImportant[i], country == country_flag_DemImportant[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(0, 1, 0.2), limits = c(-0.1,1.1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      geom_errorbar(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, ymin = LL, ymax = UL), colour = cols[1], width=.5, position=position_dodge(0.1)) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020),limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}


for(i in 1:length(predList_Essentials)){
  
  message(i)
  
  pg_list_Essentials[[i]] <- local({
    i <- i 
    pp.data <- predBornAdult_Essentials %>% filter(dv_name == name_flag_Essentials[i], country == country_flag_Essentials[i])
    pp_period.data <- predPeriod_Essentials %>% filter(dv_name == name_flag_Essentials[i], country == country_flag_Essentials[i])
    p1 <- ggplot() + 
      geom_errorbar(data = pp.data, aes(x=born_adult, ymin = LL, ymax = UL), alpha = 0.2, colour = cols[2]) +
      geom_point(data = pp.data, aes(x=born_adult, y=fit), alpha = 0.1, size = 1, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = fit), se = F, alpha = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = LL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) + 
      geom_smooth(data = pp.data, aes(x=born_adult,y = UL), se = F, alpha = 0.5, linetype = "solid", size = 0.5, colour = cols[2]) +
      ggtitle(paste(pp.data$country[1])) + 
      scale_y_continuous(name="", breaks=seq(0, 1, 0.2), limits = c(-0.1,1.1)) + 
      geom_point(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit), color = cols[1]) + 
      geom_text(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, y = fit, label=year),vjust = 1, hjust = 1, size = 2) + 
      geom_errorbar(data = pp_period.data, aes(x = (((year - 1980)*100)/40) + 1920, ymin = LL, ymax = UL), colour = cols[1], width=.5, position=position_dodge(0.1)) + 
      scale_x_continuous(name = " ",breaks = c(1920, 1940, 1960, 1980, 2000, 2020),limits = c(1920, 2020), sec.axis = sec_axis(~(((. - 1920)*40)/100) + 1980, breaks = c(1980, 2000, 2020), name = "") ) + 
      theme(
        axis.line = element_line(colour = "black"), 
        text = element_text(size=7), 
        axis.text = element_text(size = 7),
        plot.title = element_text(margin = margin(0,0,0,0), size = 7, vjust = -4),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt")
      ) + 
      theme_light() + 
      theme(
        plot.title = element_text(margin = margin(0,0,-5,0)),
        plot.margin = unit(c(0.5,5.5,0.5,5.5), "pt"),
        axis.text.x.bottom = element_text(angle = -90, hjust = 0, vjust = 0)
      )
    print(p1)
  })
  
  
}





for(i in 1){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_GAM_cohort_DemImportant.png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list_DemImportant[[i]], pg_list_DemImportant[[i+1]], pg_list_DemImportant[[i+2]], pg_list_DemImportant[[i+3]], pg_list_DemImportant[[i+4]], pg_list_DemImportant[[i+5]], pg_list_DemImportant[[i+6]], pg_list_DemImportant[[i+7]], pg_list_DemImportant[[i+8]], pg_list_DemImportant[[i+9]], pg_list_DemImportant[[i+10]], pg_list_DemImportant[[i+11]], pg_list_DemImportant[[i+12]], pg_list_DemImportant[[i+13]], pg_list_DemImportant[[i+14]], pg_list_DemImportant[[i+15]], nrow = 4))
  dev.off()
}


for(i in seq(1,length(pg_list_Essentials), by = 15)){
  agg_png(file=paste0(wdShinyPlots, "plotsSub_GAM_cohort_Essentials_",i,".png"), width = 3360, height = 2240, res = 300)
  
  grid.arrange(arrangeGrob(pg_list_Essentials[[i]], pg_list_Essentials[[i+1]], pg_list_Essentials[[i+2]], pg_list_Essentials[[i+3]], pg_list_Essentials[[i+4]], pg_list_Essentials[[i+5]], pg_list_Essentials[[i+6]], pg_list_Essentials[[i+7]], pg_list_Essentials[[i+8]], pg_list_Essentials[[i+9]], pg_list_Essentials[[i+10]], pg_list_Essentials[[i+11]], pg_list_Essentials[[i+12]], pg_list_Essentials[[i+13]], pg_list_Essentials[[i+14]], nrow = 4))
  dev.off()
}


